# include "scene.h"

#ifdef SOLVE_ENERGY

# include "discr.h"
# include <stdio.h>

#ifdef _OMP_TEMP_ADVANCE_
# include <omp.h>
#endif

# ifdef THMCOND_TDEP
extern double thm_cond(double T);
# endif

extern void temp_bc(struct grid *g, double ****T0, short tflag, unsigned short iRK);
// extern void temp_bc(const struct grid *g, double ****T, short tflag, unsigned short iRK/*, double *pE_accum*/);

/* Solving the Energy Equation using 4-stage RK-3 scheme. */
extern double temp_advance(struct grid *g, const struct arr_struct *depVar)
  {
   short i, j, k, nr, nth, nz, isec;
   double ***_T[Nsec][5], ***u, ***v, ***w, ***T, ***Tnew, ***T0, ***_T0[Nsec], *r, *th, *z, density0, E_absr;
   char vhsrc_flag;
#  ifndef  THMCOND_TDEP
   double thmcond0;
#  endif

   density0 = g->density0;
#  ifndef  THMCOND_TDEP
   thmcond0 = g->thmcond0;
#  endif

   for(isec = 0; isec < Nsec; isec++)
     {/* Assigning space for different RK-stages. */
      _T[isec][0] = g->T[isec];
      _T[isec][1] = depVar[isec].q[7]; _T[isec][3] = depVar[isec].q[7];
      _T[isec][2] = depVar[isec].q[8]; _T[isec][4] = depVar[isec].q[8];

      _T0[isec] = _T[isec][0];
     }

   nth = g->nth; th = g->th;

   if(g->vhsrc_flag == 'y')
      {
    /* If temperature in Section-1 reaches the Expected temperature rise, turn-off the volume heat source. */
       isec = 0;
       if( g->T[isec][NrG+2*g->nr[isec]/3][NthG+nth/2][NzG+g->nz[isec]/2]-T_INLET > EXPECTED_T_RISE) g->vhsrc_flag = 'n';
      }
   vhsrc_flag = g->vhsrc_flag;

/* ------------------- Calculating T(1) ------------------- */
   for(isec = 0; isec < Nsec; isec++)
       {
     /* These are used in the evaluation of E_ijk. */
        nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
        u = g->um[isec]; v = g->vm[isec]; w = g->wm[isec];
        T = _T[isec][0]; Tnew = _T[isec][1]; T0 = _T[isec][0];

#       ifdef _OMP_TEMP_ADVANCE_
        #pragma omp parallel for private(i,j,k) schedule(dynamic,1)
#       endif
        for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++)
        for(k = NzG; k < nz+NzG; k++) Tnew[i][j][k] = T[i][j][k] + 0.5*g->dt*E_ijk;

        g->T[isec] = Tnew;
       }
   temp_bc(g, _T0, 0, 1);
/* -------------------------------------------------------- */

/* ------------------- Calculating T(2) ------------------- */
   for(isec = 0; isec < Nsec; isec++)
       {
     /* These are used in the evaluation of E_ijk. */
        nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
        u = g->um[isec]; v = g->vm[isec]; w = g->wm[isec];
        T = _T[isec][1]; Tnew = _T[isec][2]; T0 = _T[isec][0];

#       ifdef _OMP_TEMP_ADVANCE_
        #pragma omp parallel for private(i,j,k) schedule(dynamic,1)
#       endif
        for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
        for(k = NzG; k < nz + NzG; k++) Tnew[i][j][k] = T[i][j][k] + 0.5*g->dt*E_ijk;

        g->T[isec] = Tnew;
       }
   temp_bc(g, _T0, 0, 2);
/* -------------------------------------------------------- */

/* ------------------- Calculating T(3) ------------------- */
   for(isec = 0; isec < Nsec; isec++)
       {
     /* These are used in the evaluation of E_ijk. */
        nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
        u = g->um[isec]; v = g->vm[isec]; w = g->wm[isec];
        T = _T[isec][2]; Tnew = _T[isec][3]; T0 = _T[isec][0];

#       ifdef _OMP_TEMP_ADVANCE_
        #pragma omp parallel for private(i,j,k) schedule(dynamic,1)
#       endif
        for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
        for(k = NzG; k < nz + NzG; k++) Tnew[i][j][k] = (2.0/3)*T0[i][j][k] + (1.0/3)*T[i][j][k] + (1.0/6)*g->dt*E_ijk;

        g->T[isec] = Tnew;
       }
   temp_bc(g, _T0, 0, 3);
/* -------------------------------------------------------- */

/* ------- Getting "T" at next time-step. It will be given by "_T[4]" --------- */
   for(isec = 0; isec < Nsec; isec++)
       {
     /* These are used in the evaluation of E_ijk. */
        nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
        u = g->um[isec]; v = g->vm[isec]; w = g->wm[isec];
        T = _T[isec][3]; Tnew = _T[isec][4]; T0 = _T[isec][0];

#       ifdef _OMP_TEMP_ADVANCE_
        #pragma omp parallel for private(i,j,k) schedule(dynamic,1)
#       endif
        for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
        for(k = NzG; k < nz + NzG; k++) Tnew[i][j][k] = T[i][j][k] + 0.5*g->dt*E_ijk;

        g->T[isec] = Tnew;
       }
   temp_bc(g, _T0, 1, 0);
/* ------------------------------------------------------------------------------ */


/* Determining total energy absorbed by the system. */
   for(E_absr = 0, isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       T = _T[isec][0]; Tnew = _T[isec][4];

#      ifdef _OMP_TEMP_ADVANCE_
       #pragma omp parallel private(i,j,k)
#      endif
          {
#          ifdef _OMP_TEMP_ADVANCE_
           #pragma omp for schedule(dynamic,3) reduction(+:E_absr)
#          endif
           for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++) for(k = NzG+nz-1; k >= NzG; k--)
           E_absr += (r[i+1]+r[i])*(r[i+1]-r[i])*(th[j+1]-th[j])*(z[k+1]-z[k])*(Tnew[i][j][k] - T[i][j][k]);

#          ifdef _OMP_TEMP_ADVANCE_
           #pragma omp for schedule(dynamic,3)
#          endif
           for(i = 0; i < nr + 2*NrG; i++) for(j = 0; j < nth + 2*NthG; j++) for(k = 0; k < nz + 2*NzG; k++)
           T[i][j][k] = Tnew[i][j][k]; /* Update the temperature with new value. */
          }

       g->T[isec] = T;
      }

   E_absr *= 0.5*g->density0*Cp_WATER;

   return E_absr;
  }
#endif
